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Abstract 



Different relaxation approximations to partial differential equations, 
including conservation laws, Hamilton- Jacobi equations, convection-diffusion 
problems, gas dynamics problems, have been recently proposed. The 
present paper focuses onto diffusive relaxed schemes for the numerical 
approximation of nonlinear reaction diffusion equations. We choose here 
a nonstandard relaxation scheme that allow the treatment of diffusion 
equations in their nonconservative form. A comparison with the tradi- 
tional approach in the case of conservative equations is also included. 
High order methods are obtained by coupling ENO and WENO schemes 
for space discretization with IMEX schemes for time integration, where 
the implicit part can be explicitly solved at a linear cost. To illustrate the 
high accuracy and good properties of the proposed numerical schemes, 
also in the degenerate case, we consider various examples in one and two 
dimensions: the Fisher- Kolmogoroff equation, the porous-Fisher equation 
and the porous medium equation with strong absorption. Moreover we 
show a test on a system of PDEs that describe an ecological model for the 
dispersal and settling of animal populations. 

1 Introduction 

The main purpose of this work is to approximate solutions of a nonlinear, pos- 
sibly degenerate, reaction-diffusion equation of the form 



for a; g 17 C R'*, d> 1, t > 0, with initial condition, u{x, 0) = uo{x) and with 
suitable boundary conditions, where the function A is a non-negative bounded 
function defined on K. The equation is degenerate if A{0) — and in this case 
the solutions often become non-smooth in finite time, developing fronts and 




(1) 



discontinuities [Aro70| . Finally, the coefRcient D is a diffusivity coefRcient and 
the function g(u) is the reaction term. 

Equations like ([T]) are relevant in describing biological and physical pro- 
cesses and are involved in the modelling of population growth and dispersal 
|GRSGM00l ITLB07a| , of waves of concentration of chemical substances in liv- 
ing organisms, of the motion of viscous fluids (see e.g. [Mur02| V Also, systems 
of equations like ([1]), coupled via the reaction terms, are capable of modelling the 
cyclic Belousov-Zhabotinskii reactions and the pattern formation on the wings 
of butterflies and the coat of mammals (see |Mur03j and references therein). 

This paper is organized as follows. In section [2 we introduce the relaxation 
approximation of nonlinear diffusion problems, in section [3] we describe the fully 
discrete relaxed numerical scheme in the reaction-diffusion case. In section |4] 
we report several numerical tests, in one and two space dimensions. The last 
application shown is about a relevant ecological model. 



2 Relaxation approximation of nonlinear diffu- 
sion 

The schemes proposed in the present work are a generalization of those proposed 
and studied in ^NP98, , and fit in the wider class of the schemes associated to 
relaxation approximations of partial differential equations, started after the well- 
known case studied in [JX95j for hyperbolic conservation laws. In the case of 
the nonlinear diffusion operator in ([T]), an additional variable v[x,t) S M'^ and 
a positive parameter e are introduced and the following relaxation system is 
considered: ^ 

{— -f div(t/) = g{u) 
(2) 
— + —A{u)Vu = V 
ot e e 

Now, formally, in the small relaxation limit, e — > 0^, system ([2]) approximates 
to leading order equation ([T]). 

We point out that in this relaxation system we do not need a third equation 
to "relax" the non-linearity of A, as was the case in ^CNPSOJ . 

In the previous system the parameter e has physical dimensions of time and 
represents the so-called relaxation time. Furthermore, each component of v has 
the dimension of u times a velocity and (/? is a velocity. The inverse of e gives 
the rate at which v decays onto —A{u)\Iu in the evolution of the variable v 
governed by the stiff second equation of ([2|). 

In order to justify the convergence of the proposed relaxation scheme to 
equation ([T]), we study the Chapman-Enskog expansion of the solutions of ([2]). 
We consider for simplicity the case with only one spatial variable and temporar- 
ily assume that A is sufficiently regular. From the second equation of ^ we 
obtain 

V = -A{u)u.^ - evt = -A{u)u.j,. - £{A{u)u.^)t + O(e^) 
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Using the relation above and ^ again, one computes 




The idea of employing a semilinear hyperbolic relaxation system to derive 
numerical schemes for (degenerate) diffusion equation was first introduced in 
|NPOO| for the purely diffusive case in the conservative form ut = A{p{u)). 
High order schemes based on such approach have been studied in |CNPS07] . As 
in that paper we introduce a parameter ip which allows to move the stiff terms 
Y^p{u) to the right hand side, without losing the hyperbolicity of the system: 



Now, the characteristic velocities of the hyperbolic left hand side are given by 
±(p and arc no longer stiff. 

We point out that degenerate parabolic equations often model physical sit- 
uations where free boundaries and discontinuities are relevant: we expect that 
schemes for hyperbolic systems will be able to reproduce faithfully these details 
of the solution. One of the main properties of ^ consists in the semilinearity 
of the system, that is all the nonlinear terms are in the (stiff) source terms, 
while the differential operator is linear. Hence, the solution of the convective 
part requires neither Riemann solvers nor the computation of the characteristic 
structure at each time step, since the eigenstructure of the system is constant 
in time. Moreover, the relaxation approximation does not exploit the form of 
the nonlinear function A and hence it gives rise to a numerical scheme that, to 
a large extent, is independent of it, resulting in a very versatile tool. 

We also anticipate here that, in the relaxed case (i.e. e = 0), the stiff source 
terms can be integrated solving a system that is already in triangular form and 
then it does not require iterative solvers. 

Obviously, when ([1]) can be put in the form ut + A(p(u)) = g{u), the re- 
laxation approximation of the parabolic part described in [CNPS07] can be 
employed. In the following sections we will compare the two possibilities from 
a theoretical and numerical point of view. 

3 Relaxed numerical schemes 

The numerical approximation of system ([3]) can be obtained following the ideas 
already exploited in |CNPS07] : first obtain a semidiscrete scheme applying an 




(3) 
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IMEX time integrator and then choose a high order non-oscillatory spatial dis- 
cretization. The reaction term will be included in the explicit part of the IMEX 
scheme. This obviously is convenient only for non-stiff reaction terms. 

For simplicity, here we describe the one-dimensional case, the generalization 
being straightforward. 



3.1 Relaxed IMEX schemes 

We observe that system ^ is in the form 

dfiz) 



dz 



da 



g{z) + -h{z) 



where z — {u,v)'^, f{z) — {v,ip'^u)'^ , g[z) — {g{u),0)'^ and h{z) = {0,—v + 
eip^Ux — DA{u)uxY' ■ When e is small, the presence of both non-stiff and stiff 
terms, suggests the use of IMEX schemes |ARS971 IKC03[ IPR05] . In the prob- 
lems considered here, the reaction term g{u) is not stiff and hence we treat it 
within the explicit portion of the IMEX scheme. 

First we semi-discretize the equation in time. Let's assume for simplicity 
a uniform time step and denote with z" the numerical approximation of 
the variable z at time f„ = nAt, for n = 0, 1, . . . We employ a j/-stages IMEX 
scheme including the 1/e terms in the implicit part, giving 



Ai^6, 



9/, 
dx 



Ai 



(4) 



where the stage values are computed as 



At 



-ai^ih{z 



(5) 



for 



fc=i 



5/, 
dx 



(fc)N 



-t- 



At 



E 

k=l 



a.^kHz^''^ 



(6) 



Here {aik,bi) and {aik,bi) are a pair of Butcher's tableaux |HNW93| of, respec- 
tively, a diagonally implicit and an explicit Runge-Kutta scheme. 

In this work we use the so-called relaxed schemes, that are obtained by 
letting e ^ in (|4]). The first stage, that is ([5]) with i = 1, is defined by 



lim ■ 











h ( 


■ yd) ■ 




{ 








—ai,ih 1 




)} 







(7) 



and thus 



,(1) 



/I) =_A(w(i)) 



dx 
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Note that, in particular, h{z^^'') = 0. Now the second stage, i — 2, reads 



.(2) 



Aia2,i 



9/ 
dx 



Hence the second component of z^^'' is determined by the stiff term of the above 
expression, namely = 0. On the other hand, due the form of h{z), there 

are no stiff terms in the equation for the first component u^'^\ which is then 
determined by a balance law. 

Summarizing, the relaxed scheme yields an alternation of relaxation steps 



h{z 







I.e. V 



-A{u 



dx 



and transport steps where we advance for time j,At 

dz df{z) 



dt 



dx 



.9(2) 



(8) 



(9) 



with initial data z ~ z^'^\ retain only the first component and assign it to u^'^'^'^\ 
Finally the value of is computed as m" + ^ biU^^\ 



3.2 Spatial reconstructions 

In order to have a fully discrete scheme, we still need to specify the space 
discretization. We will use discretizations based on finite differences, in order 
to avoid cell coupling due to the source terms. 

Recall that the IMEX technique reduces the integration to a cascade of 
relaxation and transport steps. The former are the implicit parts of ^ and ([5|), 
while the transport steps appear in the evaluation of the explicit terms i?'*^ in 
Since ([7]) and ^ involve only local operations, the main task of the space 
discretization is the evaluation of d^f, where we will exploit the linearity of / 
in its arguments. 

In the one-dimensional case, let us consider a uniform grid on [a, b] C M, 
Xj — a ~ ^ + jh for j = 1, ... ,m, where h = [h — a)/m is the grid spacing and 
m the number of cells. We denote with z" the value of the quantity z at time t" 
at Xj, the centre of the j^^ computational cell. The fully discrete scheme may 
be written as 

-^?-^t±k {f^U - ^:^r/.) + V E b^^f)^ 

1=1 1=1 

(i) 

where ^^'^ numerical fluxes, which are the only item that we still 

need to specify. It is necessary to write the scheme in conservation form and 
thus, following |S089| . we introduce the function F such that 

1 px+h/2 

f{z{x,t))^- F{s,t)As 

Jx-h/2 
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and hence 

dx 



The numerical flux function approximates F{xj^i/2)- 

In order to compute the numerical fluxes, for each stage value, we reconstruct 
boundary extrapolated data z^^-^^^ with a non-oscillatory interpolation method, 

starting from the point values ^j*-* of the variables at the centre of the cells. Next 
we apply a monotone numerical flux to these boundary extrapolated data. 

To minimize numerical viscosity we choose the Godunov flux, which, in the 
present case of a linear system of equations, reduces to the upwind flux. In 
order to select the upwind direction we write the linear system with constant 
coefficients ([9]) in characteristic form. The characteristic variables relative to 
the eigenvalues if, —if are respectively 

^ _ V + ipu ^ ^ ipu — V 
Note that u = C/ + V". Therefore the numerical flux in characteristic variables 

is + = (</5t^7+i/2,-</'V^>l/2)- 

The accuracy of the scheme depends on the accuracy of the reconstruction 
of the boundary extrapolated data. For a first order scheme we use a piecewise 
constant reconstruction such that U^j^-^j^ = Uj and Vj^ip = ^j+i- For higher 
order schemes, we use ENO or WENO reconstructions of appropriate accuracy 
[S089] . 

Since the transport steps need to be applied only to we have 
fc=l 

Finally, taking the last stage value and going back to conservative variables, 
,,"+1 - _ A h- Tf?/'^" +7/'^+ - C?/'^" ^1 



+ - % + l/2 " (^j-1/2 - ™j-l/2)] ) 



We wish to emphasize that the scheme reduces to the time advancement 
of the single variable u. Although the scheme is based on a system of three 
equations, the construction is used only to select the correct upwinding for the 
fluxes of the relaxed scheme and the computational cost of each time step is not 
affected. 



3.3 Numerical scheme 

Employing a Runge-Kutta IMEX scheme of order p, can give an integration 
procedure for ([1]) which is of order up to 2p with respect to h (see [CNPSOT] ). 
because the CFL restrictions are of parabolic type (At < Ch?). We observed 
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that this theoretical convergence rate can be achieved in practice, with careful 
choice of the approximations of the spatial derivatives in the discretization of 



Summarizing, the relaxed schemes that we propose for the numerical integra- 
tion of ((D) consists of the following steps. For each Runge-Kutta stage we need 
to compute the variables w*^*-' (relaxation steps). This requires the approxima- 
tion of a spatial gradient operator, for which we choose a central finite difference 
operator of order at least 2p. Close to the boundary of the domain, there would 
not be enough information to compute the gradient with high order centered 
finite difference formulas. We employ instead asymmetric formulas of order 2p. 
Then we need to solve the transport equation by diagonalizing the linear system 
([9]), reconstructing the characteristic variables U^^^ and y^*^ at cell boundaries 
and computing the fluxes. Again, we must choose a spatial reconstructions of 
order at least 2p and, to avoid spurious oscillations, we employ ENO or WENO 
non-oscillatory procedures. These procedures compare, for each cell, the recon- 
structions obtained using different stencils and choose the least oscillatory one 
(ENO) or compute a weighted linear combination of all of them (WENO). For 
details, see |S089) . 

Boundary conditions are enforced by extrapolating the approximate solution 
and/or the stage values u^- •* to 1 ghost point located outside the domain before 
the computation of each Runge-Kutta stage. This is done via a polynomial p{x) 
of degree 2p, fitting the first 2p values of u inside the domain and satisfying the 
boudary condition. 

4 Numerical results 

4.1 Convergence 

We tested the convergence rate of our schemes with smooth initial data. We 
considered equation ([T|) with A{u) = 1 and D = 1 and g{u) — 0. We computed 
the numerical solution with the relaxed schemes, using various combinations 
of spatial reconstructions and time integrators, starting from the initial data 
^(a:;) = sin(27ra;) for x £ [0, 1]. The errors and convergence rates are reported 
in Table [TJ This test shows that the schemes can reach the predicted order of 
convergence. 

4.2 Travelling waves 

We consider the following generalization of the Fisher-Kolmogoroff equation 



The existence of travelling waves can be proved for a wide range of param- 
eters p, q, m [Mur02j . For example in the paper [New83| gives the expression of 



m 




(10) 
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N=12 


N=36 


N=108 


N=324 


N=972 


EN02, RKl 


0.00068189 


0.00024698 


5.9722C-005 


1.1944e-005 


9.4124e-007 


EN03, RK2 


0.00067142 


0.00023703 


1.2878O-005 


8.8231e-007 


3.2784e-008 


EN04, RK2 


0.00066558 


0.00024014 


9.7088e-006 


2.2072C-007 


2.0655e-009 


EN05, RK3 


0.00066357 


0.00022958 


4.2855e-006 


3.3662C-008 


1.5442e-010 


EN06, RK3 


0.0006615 


0.0002267 


4.2338e-006 


1.0991e-008 


1.833e-011 


WEN03, RK2 


0.0006839 


0.00014311 


7.484e-005 


1.014e-005 


4.4262e-007 


WEN05, RK3 


0.00066849 


0.0001484 


2.387e-005 


2.3906e-007 


3.7223e-010 





N=36 


N=108 


N=324 


N=972 


EN02, RKl 


0.92439 


1.2922 


1.465 


2.3127 


EN03, RK2 


0.94777 


2.6512 


2.4401 


2.997 


EN04, RK2 


0.92793 


2.9202 


3.4442 


4.2522 


EN05, RK3 


0.96613 


3.6237 


4.4116 


4.9011 


EN06, RK3 


0.97477 


3.6232 


5.4194 


5.8222 


WEN03, RK2 


1.4238 


0.59004 


1.8194 


2.8504 


WEN05, RK3 


1.37 


1.6633 


4.1904 


5.8847 



Table 1: Errore in norma 1 fatte suUa soluzione con 2916 punti e velocit di 
convergenza 




Figure 1: Travelling waves with a = 2 (left) and a = 5 (right). The circles 
represent the numerical solution at t = 0, 1, 2, 3, 4, 5 and the solid lines are the 
corresponding exact solutions (fTT|) . 
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a = 2 


a = 5 




Section [3] 


CNPS07) 


Section [3] 


CNPS07, 


N 


\\Eh 


rate 


\\Eh 


rate 


l|i^l|2 


rate 


\\Eh 


rate 


30 


2.3515C-01 




6.9686e-02 




4.5251e-01 




5.0452e-02 




60 


6.9899e-02 


-1.75 


4.7352e-03 


-3.87 


2.0266e-01 


-1.15 


2.7865e-02 


-0.85 


120 


2.4679e-02 


-1.50 


4.1303e-03 


-0.19 


9.0192e-02 


-1.16 


1.3893e-02 


-1.00 


240 


6.03856-03 


-2.03 


7.0980e-04 


-2.54 


3.2529e-02 


-1.47 


3.1202e-03 


-2.15 


480 


2.1158e-03 


-1.51 


3.2753e-04 


-1.11 


1.2977e-02 


-1.32 


1.1481e-03 


-1.44 


avg. rate 




-1.71 




-1.82 




-1.28 




-1.40 



Table 2: Errors and convergence rates on the travelling wave solution pT|) . For 
a = 2 and a = 5 we compare the results of the schemes proposed in this paper 
with those described in |CNPS07j . The last row is the average convergence rate. 



such a wave for the case p = 1 and q = m = a. The speed c and the exact 
solution are: 



1 



VT+a 



i{t, x) = 



(l-e^) 



(11) 



Figure [T] shows the numerical solutions obtained for t — 5, using 300 points 
in the interval x € [—5,5]. The initial data and the exact solutions shown are 
defined according to (fTTj) . The accuracy of the schemes is good in both cases. 
In the case with a = 5 there is a lower accuracy around the corner, due to the 
higher degeneracy of the diffusion term and the more pronounced stiffness of 
the reaction term. 

Table [2] shows the errors and convergence rates, measured at t — 5, compar- 
ing the results of the schemes presented in this paper and those obtained with 
the schemes described in [CNPS07' . In order to apply the schemes of |CNPS07] , 
equation pop was rewritten in the "conservative form" ut = [u"^ /m]xx + ^(1 — 
It"), where m ~ a + 1. We point out that, although we employed the Heun 
scheme and ENO reconstruztions of order 3, the convergence rate is limited by 
the C*^°-' regularity of the solution. The comparison shows that the schemes 
for the conservative form perform slightly better than those described here for 
the nonconservative form, both regarding the convergence rate and the absolute 
value of the errors. Both schemes could be improved by integrating the reaction 
term with the diagonnally implicit part of the IMEX scheme, expecially in the 
test with a = 5. 

Thus, whenever the equation can be put in the conservative form required 
by the schemes of [CNPSOT] . those are to be preferred. In the other cases, 
which frequently occour in more than one dimensions or in the case of systems 
of equations, the present scheme is a good choice for the numerical integration 
of the differential problem. 
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4.3 Two dimensional tests 



We tested the multi-dimensional version of our integration scheme on the equa- 
tion 

— ^A{un-cuP (12) 

for X € [—2, 2]^ and t > 0. The above equation presents interesting finite-time 
extinction phenomena, reported in |Mik95j . 

We tested the ability to reproduce the extinction phenomenon and the persis- 
tence of asymmetry in the initial datum along the evolution. We took u{x, y, 0) 
as a radially symmetric function with a small perturbation, see Figure [5J and 
evolved it with the same equation and parameters as before, until extinction. 
Figure [2] shows clearly that the initial perturbation of the radial symmetry is 
maintained until the solution vanishes. 



4.4 Application to a system of PDEs 

We consider now an ecological model described in |TLB07b] . It consists in a 
system of PDEs describing the dispersal and settling of frogs. The model splits 
the population of frogs in two classes: the ones that have already settled and 
those that are still migrating. The moving class diffuses with rate D and settles 
at rate S. Once an individual is in the settled class, they cannot migrate any 
more. The authors of [TLB07b] study various mechanisms of dispersal and 
settling: we describe here only those that they think relevant in the case study 
considered in the paper, that is density-dependent dispersal and constant or 
u-dependent settling rate. 

Let Um (t, x) and Us (t, x) be the population density of migrating and settled 
animals, D{u) and S{u) suitable functions describing the dispersal and settling 
mechanism. The population is described by the following parabolic equations. 

yS/.^ ■ {D{Ujn + Us)Va;(Wm)) ~ S{Ujn + Us)Um (13a) 
S{Ujn+Us)Um (13b) 

Initial conditions mimic a population concentrated around x = for Um(0,x) 
and Us(0,x) = 0. In the case D(u) — u/uq (for constant uq) considered in the 
paper, the diffusion equation (|13a[) is degenerate and gives rise to Barenblatt- 
like profiles for the density Um of the moving population. We point out that 
a very accurate numerical approximation is very important, since there is no 
mechanism in the model to move the population once it has settled. 

If a second population v is released in the same spot, the model is augmented 



dum 
dt 

dus 

'dt 
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t=0 50046 



t=5.0007 





Figure 3: Model ^ for frog populations in [TLB07b]. Solutions at t = 0.5, 
t = 5, t = 5.5, t = 10, t = 15 and t — 20. Note the different scale in the top 
right panel. In blue: the u population released at x — at t = 0; in green: the 
V population, that is released at a; = after the first has settled {t — 5 here). 
The black dashed line is the total settled population, i.e. Ug + Vg- 
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with the following equations: 

aAx{cu) + f3{um + Us - Cu) (13c) 



fiVx ■ {D{Vm + Vs)Vx{Vm)) + x ' [Vm'^ x{Cit)) - S{Us + + Vs)Vm 

(13d) 



dcu 
dt 

dt 

dv 

^ S{Us+V„i+Vs)Ujn (13e) 

Here Vm and Vs are the densities of migrating and settled animals belonging to 
the second population, while c„ is a diffusible pheromone released by the first 
population and discouraging the v population from staying in the area where 
the u density is higher. 

The case study of |TLB07b] is about a release of a batch of frogs (v) at 
the same spot where an earlier release (u) had taken place, at a time when the 
frogs from the first batch have already settled down. In order to reproduce the 
situation of the case study of [TLB07b] . one has to evolve in time the system 
(jl3p starting with an initial nonzero Um population while Us = Vm = fs = and 
adding a nonzero Vm density after the u population had time to settle settled 
(that is when Um is negligible). 

We employed free-flow boundary conditions and chose parameters as sug- 
gested in ITLlBOTa) : 

fi^l 7=^0 a = 0.01 /3==10 

The initial population is taken to be 

and we chose 

^(")-ol5 ^("'^^-^(l-ols) 
where xi^) — 1 when x > and xi^) = otherwise. 

The u population settles in the time lapse from t = to i = 5 with maximum 
density located some distance apart from the original release spot (x = in the 
simulation). After this, the release of the second batch is simulated setting, at 
t — 5, Wm(5, x) — Um{0, x). Thus the same initial condition is employed for both 
batches. 

Figure [3] clearly shows that this second population first moves further away 
than the peak of the Ug distribution and start settling there. Only afterwards 
the V population settles also near x = 0. 

The panel in figure [3] shows the evolution of the system, from the early u 
dispersal that gives rise to a Barenblatt-like profile for Um (upper left), to the 
release of the v population when Um — (upper right), the early migration and 
settling of the second population around x = 2 (middle), the beginning of the 
secondary settling at a; = (bottom left) and the final steady state (bottom 
right). 
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We reproduced with our method also all the other cases studied in |TLB07aj . 
We point out that the standard integration procedure used for the simulations 
in |TLB07a] required a smoothing of the x function and, despite that, produces 
unrealistic wiggles in the population densities, that are not seen with the method 
presented here. 



5 Conclusions 

We have proposed and analyzed relaxed schemes for nonlinear degenerate re- 
action diffusion equations. By using suitable discretization in space and time, 
namely ENO/WENO non-oscillatory reconstructions for numerical fluxes and 
IMEX Runge-Kutta schemes for time integration, we have obtained a class of 
high order schemes. The theoretical convergence analysis for the semidiscrete 
scheme and the stability for the fully discrete schemes have been studied by us, 
for the case of nonlinear diffusion, in [CNPS07] . 

Here we tested these schemes on travelling waves solutions and in cases 
where the solution vanishes in finite time. In all cases we observed a very good 
agreement with known properties of the exact solutions. 

Moreover we have shown how the accurate approximation of the fronts aris- 
ing in degenerate parabolic equations can lead to improved accuracy in the 
simulation of an ecologically relevant model. 
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